# annual yearly anaysis 
# hypothesis: avg temp VS anchovy catch
# winter temp vs anchovy catch
# try ccf, offset, etc
library("ncdf4")
library("raster")
library("tidyverse")
library("readxl")


sst <- stack("input/cmems_mod_blk_phy-tem_my_2.5km_P1M-m_1715770301759.nc")
anchovy <- read_excel("input/bs_anchovy.xlsx")

years <- 1993:2021
avgs <- c()
sd <- c()
for (y in years) {
  idx <- y - 1993
  
  month.st <- idx * 12 +1
  month.end <- idx*12 +12
  #print(paste0("st= ", month.st, ", end=", month.end))
  
  sst.slice <- sst[[(month.st:month.end)]]
  sst.slmean <- mean(sst.slice)

  avgs <- c(avgs, as.numeric(cellStats(sst.slmean, "mean")))
  sd <- c(sd, as.numeric(cellStats(sst.slmean, "sd")))
}

dfan <- data.frame(year = years, mean = avgs, sd = sd)
for (i in 1:nrow(dfan)) {
  row <- dfan[i,]
  catch <- anchovy[which(anchovy$year == row$year), "crimea_ua_nw"]
  dfan[i, "catch"] <- as.numeric(catch)
}

png("output/plots/ccf-annualsst-catch.png", res=180, width=1200, height=900)
ccf(dfan$mean, dfan$catch, main = "Корреляция Спирмана среднегодовые ТПВ и уловы хамсы", lag.max=5)
dev.off()

anualcf <- ccf(dfan$mean, dfan$catch, lag.max = 5)
write.csv(anualcf$acf, "output/anual-ccf.csv", row.names = F)

years <- 1993:2021
avgs <- c()
sd <- c()
for (y in years) {
  idx <- y - 1993
  
  month.st <- idx * 12 +12
  month.end <- idx*12 +14
  #print(paste0("st= ", month.st, ", end=", month.end))
  sst.slice <- sst[[(month.st:month.end)]]
  sst.slmean <- mean(sst.slice)

  avgs <- c(avgs, as.numeric(cellStats(sst.slmean, "mean")))
  sd <- c(sd, as.numeric(cellStats(sst.slmean, "sd")))
}

dfwn <- data.frame(year = years[-1], mean = avgs[-1], sd = sd[-1])
for (i in 1:nrow(dfwn)) {
  row <- dfwn[i,]
  catch <- anchovy[which(anchovy$year == row$year), "crimea_ua_nw"]
  dfwn[i, "catch"] <- as.numeric(catch)
}

png("output/plots/ccf-wintersst-catch.png", res=180, width=1200, height=900)
ccf(dfwn$mean, dfwn$catch, main = "Корреляция Спирмана среднезимние ТПВ и уловы хамсы", lag.max=5)
dev.off()

write.csv(dfwn, "output/winter-sst-catch.csv", row.names = F)

wintercf <- ccf(dfwn$mean, dfwn$catch, lag.max = 5)
write.csv(wintercf$acf, "output/winter-ccf.csv", row.names = F)
